library(dplyr)
library(ggplot2)
library(qqman)
library(ggrepel)
library(plotly)
library(manhattanly)
mhp_int <- function(data_loc) {
dat <- read.table(data_loc, header = T) %>%
filter(!is.na(TWAS.P) & TWAS.P != 0) %>%
rename(P = TWAS.P, SNP = BEST.GWAS.ID, GENE = ID) %>%
mutate(BP = (P0 + P1)/2,
LOGP = -log10(P)) %>%
select(SNP, GENE, CHR, BP, P, LOGP) %>%
group_by(CHR)
manhattanly(dat, snp = "SNP", gene = "GENE",
annotation1 = "P", annotation2 = "LOGP")
}
mhp_info <- function(data_loc, logp_bar = 3) {
dat <- read.table(data_loc, header = T) %>%
filter(!is.na(TWAS.P) & TWAS.P != 0) %>%
rename(P = TWAS.P, SNP = ID) %>%
mutate(BP = (P0 + P1)/2) %>%
select(SNP, CHR, BP, P) %>%
group_by(CHR)
dat_plot <- dat %>%
summarise(chr_len = max(BP)) %>%
mutate(tot = cumsum(chr_len) - chr_len) %>%
select(-chr_len) %>%
left_join(dat, ., by=c("CHR"="CHR")) %>%
arrange(CHR, BP) %>%
mutate(BPcum = BP + tot) %>%
mutate(is_annotate = ifelse(-log10(P) > logp_bar, "yes", "no"))
axisdf <- dat_plot %>%
group_by(CHR) %>%
summarize(center = (max(BPcum) + min(BPcum)) / 2)
ggplot(dat_plot, aes(x=BPcum, y=-log10(P))) +
geom_point(aes(color=as.factor(CHR)), alpha=0.8, size=1.3) +
scale_color_manual(values = rep(c("grey", "skyblue"), 22 )) +
scale_x_continuous(label = axisdf$CHR, breaks= axisdf$center ) +
scale_y_continuous(expand = c(0, 0) ) +
geom_label_repel(data=subset(dat_plot, is_annotate=="yes"), aes(label=SNP), size=2) +
theme_bw() +
theme(
legend.position="none",
panel.border = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()
)
}
ad_Brain_Amygdala
mhp_info("./Data/ad_Brain_Amygdala.dat")
mhp_int("./Data/ad_Brain_Amygdala.dat")
ad_Brain_Cerebellum
mhp_info("./Data/ad_Brain_Cerebellum.dat")
mhp_int("./Data/ad_Brain_Cerebellum.dat")
ad_Brain_Cortex
mhp_info("./Data/ad_Brain_Cortex.dat")
mhp_int("./Data/ad_Brain_Cortex.dat")
ad_Brain_Hippocampus
mhp_info("./Data/ad_Brain_Hippocampus.dat")
mhp_int("./Data/ad_Brain_Hippocampus.dat")
t2d_Adipose_Subcutaneous
mhp_info("./Data/t2d_Adipose_Subcutaneous.dat")
mhp_int("./Data/t2d_Adipose_Subcutaneous.dat")
t2d_Adipose_Visceral_Omentum
mhp_info("./Data/t2d_Adipose_Visceral_Omentum.dat")
mhp_int("./Data/t2d_Adipose_Visceral_Omentum.dat")
t2d_Pancreas
mhp_info("./Data/t2d_Pancreas.dat")
